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Abstract. We discuss a projector Monte Carlo method for quantum spin models 
formulated in the valence bond basis, using the S — 1/2 Heisenberg antiferromagnet 
as an example. Its singlet ground state can be projected out of an arbitrary basis 
state as the trial state, but a more rapid convergence can be obtained using a 
good variational state. As an alternative to first carrying out a time consuming 
variational Monte Carlo calculation, we show that a very good trial state can be 
generated in an iterative fashion in the course of the simulation itself. We also show 
how the properties of the valence bond basis enable calculations of quantities that 
are difficult to obtain with the standard basis of S z eigenstates. In particular, we 
discuss quantities involving finite-momentum states in the triplet sector, such as 
the dispersion relation and the spectral weight of the lowest triplet. 



1 Introduction 

Quantum Monte Carlo (QMC) simulations of spin systems have traditionally 
been carried out in the basis of eigenstates of the spin-z operators <Sf, i = 
1,...,N, i.e., the basis of "up" and "down" spins in the case of S = 1/2 
(which is the case we consider here). For the prototypical model of interacting 
quantum spins, the antifcrromagnetic (J > 0) Heisenberg hamiltonian, 

h = jJ2 ^ • a, = jJ2lsiSj + h(s+sr + srsf)], (1) 

(*>i) 

this basis is clearly natural and convenient, as an off-diagonal operator acting 
on a basis state just flips two spins or destroys the state. Starting with the 
work of Suzuki [1], finite-temperature simulation methods employing the spin- 
z basis were developed in which a quantum mechanical expectation value for 
a system in D dimensions is mapped onto an anisotropic classical statistical- 
mechanics problem in D + 1 dimensions — the discretized [2-4] or continuous 
[5,6] imaginary-time path integral. There are now very efficient methods uti- 
lizing loop-cluster [7,8,5] or "worm" [6] updates of the world-line spin config- 
urations. These methods have enabled studies of systems with w I0 4 — 10 5 
spins in the low-temperature (ground-state) limit and much more at elevated 
temperatures. Loop updates have been developed [9,10] also for the alterna- 
tive and now frequently used power-series expansion representation [11-14] 
of the partition function (stochastic series expansion; SSE), where the spin- z 
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Fig. 1. Two valence bond states \Vk), |V;) in two dimensions and their overlap graph 
corresponding to (Vj|V]fc) = 2 N °~ N ^ 2 , where iVo is the number of loops formed (in 
this case N a = 3 and the number of sites N = 16). Filled and open circles correspond 
to sublattices A and B. The sign convention in Eq (3) for a singlet valence bond 
dictates that spins i and j belong to A and B, respectively 



basis is also normally used. It is in principle possible to adapt these ap- 
proaches to other local bases, e.g., that of singlet and triplet states of spin 
pairs on a dimerized lattice. This basis is often used in diagrammatic and 
series-expansion calculations [15], but its implementation in QMC simula- 
tions is typically rather cumbersome. 

Zero-temperature (T = 0) simulations, in which the ground state is pro- 
jected out of a trial wave function, are also normally carried out in the 
spin-z basis [16,17]. Here we will discuss an alternative approach to ground- 
state calculations which turns out to have some unique features enabling 
access to quantities that are normally difficult to obtain with standard finite- 
temperature or projector methods. We make use of the valence bond basis, 
i.e., states in which the spins are paired up into singlets; 

\V) = \{h,h){i2,32) ■ ■ ■ {iN/ijN/i})- (2) 
Here (i,j) denotes a singlet formed by the spins at sites i and j; 

= (I T4i) - U<Ti»/V2, (3) 

and the total number of sites N is assumed to be even. While in principle 
one can include all possible pairings of the spins, it is in most cases better to 
consider a smaller basis in which the sites are first divided into two groups, 
A and B, of N/2 spins each, and to only consider singlets in which 

the first index i £ A and the second j £ B [18-20]. In the case of a bipartite 
lattice, these groups are naturally the two sublattices, as shown in Fig. 1. This 
restricted VB basis has (N/2)\ states and is still massively overcomplete — the 
singlet space has N\/[(N/2)[\ 2 (N/2 + 1) dimensions [18]. The VB basis states 
are all non-orthogonal, overlapping with each other according to the simple 
loop rule illustrated in Fig. 1. 

The VB basis was introduced already in the early 1930s [18,21,22] and 
has played an important role in exactly solvable models [18,23-25]. Later, 
it became a tool for describing spin liquids — the resonating valence bond 
(RVB) mechanism introduced by Fezekas and Anderson [26-28], in which the 
ground state is dominated by short valence bonds. In exact diagonalization 
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studies, the VB basis is useful in cases where it is a good approximation to 
only consider a restricted (and incomplete) space of short bonds (spin liquids 
and other states with no magnetic long-range order) [29-32] . Variational cal- 
culations in the VB basis have been carried out for the 2D Heisenberg model 
[20,33]. Furthermore, Liang realized that a variational VB state could be con- 
siderably improved by stochastically projecting it with an operator (—H) m 
for large m [34]. Later, Santoro et al. devised a Green's function method for 
calculating energies in the VB basis [35]. Despite the promising results ob- 
tained in these studies, there was, to our knowledge, no further developments 
of QMC methods in the VB basis until one of us recently introduced two 
related projector algorithms [36], improving on the schemes of Liang [34] and 
Santoro et al. [35]. These algorithms have already been applied in studies of 
quantum phase transitions [37,38] and entanglement entropy [39]. 

Some previously unnoticed advantages of the VB basis in QMC algorithms 
were pointed out in Ref. [36]. Here we summarize our recent work on VB 
projector methods and highlight some of their unique features. We discuss 
in particular a scheme for "self-optimizing" the trial state out of which the 
ground state is projected, and also show how to study properties of triplet 
excitations at finite momentum. 



2 Ground state projection 

Consider a state {iP} and its expansion in terms of eigenstates \ri), n = 0, 1, . . ., 
of some hamiltonian H; 

|#> =$>!»>■ (4) 

n 

With C a constant chosen such that the lowest eigenvalue Eq — C is the 
largest in magnitude, a large number m of repeated operations with C — H 
projects out the ground state, 



(C - H) m \V) ^ c Q (C - E y 



(5) 



provided that the overlap Co ^ 0. Here we will first consider singlet eigenstates 
of the Heisenberg model (1), which can be expanded in VB states; 

l#>=£/iM- (6) 

i 

Because of the overcompleteness of the VB basis, this expansion is not unique. 
That, however, does not prohibit the ground state |0) to be projected out 
according to Eq. (5). The Heisenberg hamiltonian can be written in terms 
of singlet projection operators Hi, = H^j/^ on the interacting spin pairs, 
{i(b),j(b)}, b = 1, . . . , Nf„ = DN (for a periodic cubic D-dimensional lattice); 

N b N b 

h = - j y~] H b = -jy~] H i(b),j(b), = -(i - Si • Sj). (7) 

b=l b=l 
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When a singlet projector Hij acts on a VB basis state, one of two things 
can happen; 1) if i,j belong to the same bond the state is unchanged with a 
matrix element 1, or 2) if they belong to different bonds these two bonds are 
reconfigured ("flipped") with matrix element 1/2; 

H i3 \---{i, ])■■■) = | ■■■(*, j)"-), (8) 
H ij \---{i,k)---(l,j)---) = l\---{i,j)---(l,k)---). (9) 

Here positive-definitness of (9) is directly related to the two sites i and j being 
in different sublattices. For a frustrated system, where there are operators 
with both sites i,j in the same sublattice, positive-definitness does not hold 
[19,20]. For a non-frustrated system the simple bond flip (9) makes for a 
convenient stochastic implementation of the ground state projection (5). We 
write the projection operator as (with J = 1 henceforth) 

(<7 - ff) m = £ #6 = £ Pr, Pr = ■ ■ ■ H b rH b r, (10) 

\6=1 / r 

where we have introduced a compact notation P, r , r = 1, . . . , N™, for the 
different strings of singlet projectors. When a string P r acts on a VB basis 
state \V) the result is another basis state, which we denote |V(r)), with 
a prefactor (weight) W r which is simply given by the number m s of off- 
diagonal operations in the course of evolving \ V) to |V(r)); 

P r \V) = W r \V(r)), W r = 2~ m °« . (11) 

We here first consider projecting the ground state out of a single VB 
basis state; later we will consider the use of a more complicated trial state. 
We consider two ways to calculate expectation values: 

(ff) Eg£l - Z r Wr(#\H\V(r)) 

i a\ _ Y.ri(v\p:AP r \v) _ Ew^r^(noi^ino) mi 

[ 1 Y,rl<y\P!Pr\V) "~ ErlWrWtiVmir)) ' 1 ' 

We will discuss how to estimate these using importance sampling; terms 
(configurations) of the numerators are illustrated in Fig. 2. We will refer to 
(12) and (13) as the single and double projection, respectively. 

In (12), which is an exact (when m — > oo) expectation value only of the 
hamiltonian (or other operators for which the ground state is an eigenstate) 
the state |<^) is in principle arbitrary. It is very convenient to use a state 
which has equal overlaps with all the VB basis states, e.g., the Neel state 
|#jv) (ah spins up on sublattice A and down on B). It is easy to see that 
(& N \V) = 2- N ' 2 for any basis state |V). Since H\V(r)) = - £ b H b \V(r)) is 
a sum of basis states multiplied by factors —1 or —1/2, the overlaps with (<^jv| 
drop out altogether and do not have to be considered further. If the projector 
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Fig. 2. Propagation of a VB state on a 6-site chain. The horizontal bars represent 
nearest-neighbor Heisenberg interactions (singlet projectors). In the single projec- 
tion (a) the state is propagated from right to left, and an estimator for the ground 
state energy is obtained by acting once more with all terms of the hamiltonian. In 
the double projection (b) the state is projected from the right and the left, and 
any operator expectation value can be estimated by calculating the corresponding 
matrix elements between the propagated states 



strings P r in (12) are importance-sampled according to their weights W r , the 
estimator for the ground state energy is thus 

E = (H) = -<m d + ±m ), (14) 

where md and m D are, respectively, the number of diagonal and off-diagonal 
operations Hb\V(r)) (and + m = Nf,). It should be noted that although 
this estimator is exact in the limit m — ► oo, it is not variational. The correct 
energy may thus be approached with increasing m from above or below. 

Eq. (13) is valid for any expectation value and in the case of A = H gives 
a variational estimate of the energy. Using W r Wi(V (l)\V (r)) as the sampling 
weight, the estimator for any (A) is of the form 

{A) = /{vmm\. (15) 



\ (V(l)\V(r)) 

In the case of a spin correlation function (Si-Sj), the matrix element is related 
to the loop structure of the overlap graph [19,20] (illustrated in Fig. 1): 



(VHS S-IV) ("+3/4' if i, j S same loop, same sublattice, 

' = < —3/4, if i,j e same loop, different sublattices, (16) 
[ 0, if i,j € different loops. 



(Vi\V r , 



Measuring the spin correlations is hence straight-forward once the overlap- 
loops have been constructed. Higher-order functions, e.g., dimer-dimer cor- 
relations ((Si • Sj)(Sfc • S;)), are also related to the loop structure [40]. 

Note again that no bond operator Hi, can destroy a VB state and that 
all the states have non-zero overlap with each other. Thus all terms in (13) 
contribute to the expectation value. This turns out to be an advantage in 
constructing a Monte Carlo algorithm, as any change made in the operator 
strings can be accepted with some probability. With an orthogonal basis, such 
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as the S? eigenstates, there would be considerable constraints, both in terms 
of individual operators in the projection [the spin flip operator in (1) can act, 
without destroying the state, only on anti-parallel spins], and in ensuring a 
non-zero overlap between the propagated states (the two propagated states 
have to be identical). Note also that the singlet projectors are non-hcrmitian 
in the VB basis. As indicated in (13), and illustrated in Fig. 2, we here 
propagate two states, \V(l)} oc Pi\V) and |V(r)) oc P r \V), and subsequently 
compute their overlap and various matrix elements. Propagating \V) with 
Pi P r and then taking the overlap with \V) is not equivalent term- by-term. 

To carry out the projection stochastically, the operator strings are stored 
in arrays [P a ] = [bf}[b%} ■ ■ • [6° ], where bf G {1, . . . , N b } with a = 1,2 cor- 
responding to P r and Pi, respectively, in the double projection; in the single 
projection a is redundant. A table holds the site pairs i(b),j(b). The state 
\V) is stored in a list [V] = [i>i][i>2] • • • [vn] where Vi = j and Vj = i if there is 
a valence bond at Propagation with the bond flips (9) is easily carried 

out in this representation. The state list [V] is then first copied into two lists, 
[Vi] and [V2], in which |V(r)) and |V(Z)) are constructed. 

The simulation can be started with a randomly generated operator string. 
The strings can can be updated in an trivial way, by changing a number R 
of operators at random. In either [Pi] or [P2], R positions pj, j = 1,...,R, 
Pj G {1, . . . ,m} (all different) are generated. Their contents 6™ are picked 
randomly from the set {1, . . . , Af,} (excluding the old value for each pj). To 
calculate the Metropolis acceptance probability, the state is propagated with 
the updated operator string and the number of off-diagonal operations m s in 
(11) is counted. In the single expansion, the acceptance probability is simply 

= min[2 m o«- m of",l], (17) 

whereas in the double projection an overlap ratio appears as well. In the 
double projection, we change operators only in one of the operator strings at 
a time, so that only one state has to be propagated. It is customary to define 
a size-normalized Monte Carlo step (or "sweep"). For projector length m we 
do m replacement attempts, and so our sweep is independent of the number 
of operators R replaced in each update. 

Normally, in Monte Carlo simulations one does not compute the full 
weight, because it is possible to more speedily calculate just the change in 
the weight [the weight ratio in (17)]. In the present formulation of the VB 
projector algorithm the weight is, however, recomputed from scratch each 
time, because a better scheme is not yet known. Each update hence requires 
on the order of m operations. In the double expansion, construction of the 
loops needed to compute the overlap scales as N, but typically m > N and 
the propagation of the state dominates the simulation. In spite of the need 
to recalculate the weight, the scheme is sufficiently efficient to compete with 
other ground state QMC methods. More importantly, as we will discuss in 
Sec. 4, the VB basis offers access to quantities out of reach for other methods. 
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Fig. 3. (a) Acceptance rate in a double projection and (b) the number of bonds 
changed in the propagated state in accepted updates for a 16 x 16 lattice (symbols 
with lines; the bare lines are for L — 8) versus the projection length (normalized 
by the number of sites N). The trial state \V) was a columnar dimer state 



The optimum number R of operators to replace depends on the accep- 
tance rate. In Fig. 3(a) we show the acceptance rate for a double-projected 
2D system versus the length m of the projection for R — 1 — 6. As expected, 
the acceptance rate decreases with increasing i?, but it does not change ap- 
preciably with m. It also depends only weakly on the lattice size. Multiplying 
the acceptance rate by R gives the average number of operators changed; it 
initially increases with R but has a maximum for R = 7 (for L = 16). The 
optimum R is clearly model/lattice dependent. Another characteristic of the 
update is the number of bonds changed in the projected state |V(r)) as a 
consequence of the modifications of P r . This number is shown in Fig. 3(b). 
It is seen to increase with R, as expected. As a function of m the num- 
ber of changed bonds decreases. This behavior reflects a loop structure of 
the singlet-projection operators [41,8], which implies that some changes "up- 
stream" in the operator string may be healed further downstream in the 
propagation. For a finite lattice in the limit m — > oo one would expect a 
substitution of an operator far upstream in P r to have no effect on the final 
propagated state |y(r)). This does not imply that this update is inconsequen- 
tial, as the sampling is over paths, not just the final states in the propagations. 
In principle, it should be possible to take advantage of the underlying loop 
structure of the singlet projectors [41,8] to devise a loop update in the VB 
basis, analogous to such updates in world-line [7,8] and SSE [9] simulations. 
However, we have not yet been able to construct a scheme that is in practice 
faster than the trivial random substitution with full state propagation. 

3 Self-optimized trial state 

So far, we have projected the ground state out of a single VB basis state \ V). 
This works well [36] , but the rate of convergence with m of course depends 
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on the state chosen; ideally one would like to maximize the overlap (V|0). 
One way to obtain a typically good single-configuration trial state is to first 
start with an arbitrary one; a regular bond pattern or a randomly generated 
configuration. After carrying out some projection steps with this state, the 
current propagated state |V^(r)) is chosen as a new trial state. Since this state 
has been generated in the projection it should contribute substantially to the 
ground state and hence typically will be better than a completely arbitrary 
one. However, as we will discuss next, we can do much better than this. 

Liang's original motivation for introducing a projector technique in the 
VB basis was to improve on a variational calculation [34] . Previously, Liang, 
Doucot, and Anderson had studied a variational amplitude-product state for 
the 2D Hcisenbcrg model of the form [20] 

N/2 

l^)=^MVfc), fk = Y[h(x bk ,y bk ), (18) 

k 6=1 

where Xbk and y b k are the x- and y-lengths of bond b in VB state k, as il- 
lustrated in Fig. 4(a). Liang et al. tried power-law and exponential forms 
depending only on the total length I of the bonds [the " Manhattan" length 
I = x + y was used, but defining I = (x 2 + y 2 ) 1 ^ 2 should not change things 
qualitatively], in addition to keeping several short-bond amplitudes as pa- 
rameters to optimize [20]. More recently, all the amplitudes were optimized 
without any assumed form on lattices with up to 32 x 32 sites, with the result 
that h(l) <~ Z~ 3 for long bonds [33]. One of us has also recently showed the 
more gener al result h(l) ~ l~( D + 1 ) within a mean-field approach for a D- 
dimensional cubic lattice [42]. In 2D the fully-optimized amplitude-product 
state turns out to be extremely good, with an energy deviating by only 
w 0.06% from the exact ground state energy and with the long-distance 
spin-spin correlations reproduced to within 2% [33]. 
With the state (18) an expectation value is given by 

/ a \ _ ^ fc P fkfp(V P \Vk} (V ty$!) 

{) ~ E rl fkMv P \v k ) ' U9j 

which can be evaluated using importance sampling of the VB configurations 
with weight f k f P (V p \Vk}- Liang et al. introduced a very simple scheme for 
updating the dimer configurations [20], which we here illustrate in Fig. 4(b). 
Choosing two next-nearest- neighbor sites, i.e., ones on a diagonal of a 4-sitc 
plaquette (or, in principle, any two sites on the same sublattice), the two 
bonds connected to them are reconfigured in the only possible way which 
maintains only bonds between the A and B sublattices, as shown in the figure. 
Labeling the two initially chosen sites 1 and 2, and the bonds connected to 
them b = 1,2, the Metropolis acceptance probability is, assuming that the 
bond update was made in |T4), resulting in \Vk>), 



-^accept 



h(xik> , yik')h(x 2 k> , V2k' ) ( V P \ Vk> ) l 
h(xik,yik)h(x2k,y2k) (V p \Vk) ' 



(20) 
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Fig. 4. (a) Definition of the size of a bond b. (b) Reconfiguration of two bonds in 
an update of the trial state. The two sites marked 1, 2 are chosen at random among 
all pairs of next-nearest neighbor sites. 



We can use an amplitude-product state as the trial state in the projector 
QMC method, using some set of amplitudes not necessarily originating from 
a variational calculation. In updating the bond configurations, we then also 
must compute the new weight of the propagation; P r \Vk>) = Wk> r \Vk> (r)), and 
of course (VJ,|Vfc) in (20) is replaced by (V p (l)\ Vfc(r)) (in the case of the double 
projection; for the single projection there is no overlap). The acceptance rate 
of a state update is similar to that of an operator string update with a small R, 
and often we find it advantageous to combine state and operator updates. To 
save some time, one can tentatively accept/reject a bond update based solely 
on an amplitude ratio — Eq. (20) without the overlaps — and then calculate 
the overlap and the propagation weight for a final accept/reject probability 
only for tentatively accepted bond updates. 

In the variational calculation (H) is minimized with respect to all h{x, y). 
With a recently developed stochastic optimization method [33] , all the cx N 
amplitudes can be minimized for moderate-size lattices (up to 32 x 32 sites 
were considered in Ref. [33]). In principle we could follow Liang [34] and use 
the best possible variational state as our trial state in the projector QMC 
method — indeed this can be expected to be the optimum starting point. 
However, we will now describe a scheme which delivers a trial state nearly as 
good as the best variational state, at a smaller computational cost. 

Consider the probability distribution P(x,y) of valence bonds. In an 
amplitude-product state (18) we would have P(x,y) ~ h(x,y), were it not 
for the "hard-core" constraint of only one bond per site. Even with this con- 
straint, it is clear that the probabilities and amplitudes are related in a mono- 
tonic way; increasing h(x,y) for some given (x,y), while keeping the other 
amplitudes fixed, will lead to a larger P(x,y). This fact can be exploited 
in constructing a good trial state. We define two different probability dis- 
tributions, Po(x,y) and P m (x,y), the former being the just discussed bond 
probability in a trial state of the form (18) and the latter the probability 
distribution in the projected state. For sufficiently large m, P m is an exact 
property of the ground state, whereas Pq is a property of the trial state and is 
in general different from P m . However, for given m, we can adjust the ampli- 
tudes h(x, y) of the trial state such that P m (x, y) — Pq(x, y) for all x, y. If this 
is done for m sufficiently large, then our trial state has a bond distribution 
identical to that of the exact ground state. Such a state is often almost as 
good as the best variational state. The reason that this is useful in practice 
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Fig. 5. Left: Energy versus projection length for a 16 x 16 system, using three 
different trial states (labeled in the inset). Right: The energy of the self-optimized 
trial state compared with the projected energy using that trial state. The dashed 
line shows the energy evaluated independently using the SSE method [43] 



is that it is very easy to adjust the amplitudes to achieve self-consistency. 
Because of the monotonous relationship between h{x,y) and Po(x, u), we can 
simply increase h(x, y) by some amount if Pq(x, y) < P m (%, y) and decrease it 
if P (x,y) > P m (x,y), and repeat this until self-consistency is achieved. We 
use the following scheme to update the amplitudes after the kth iteration; 

hx[h(x, y)} -» ln[h{x, y)] + RAN • 0{k) • sign[P m (x, y) - P Q (x, y)], (21) 

where RAN is a random number in the range [0, 1) and 0(k) decreases with 
the iteration step k = 1, 2, . . ., according to 0{k) oc k~ a . For the exponent, 
we typically use a — 3/4. To evaluate the probabilities Pq and P m (in two 
independent simulations; with and without projection of the trial state), the 
number of Monte Carlo sweeps does not have to be very large, because we 
only need the sign of the difference of the two probabilities. We normally 
use on the order of 100-1000 sweeps per iteration. Even if the stochastically 
evaluated sign in (21) occasionally may be wrong, it is correct on average 
and the amplitudes typically converge to a self-consistent solution after a 
few hundred iterations. Due to the stochastic nature of the procedure, self- 
consistency of course obtains only to within some statistical error, which can 
be reduced by increasing the number iterations and/or sweeps per iteration. 

Fig. 5 shows results for the energy of a 16 x 16 system obtained in double 
projections with three different trial states; a columnar dimer state, a ran- 
domly picked state generated while projecting the columnar state, as well as 
the self-optimized state. Already for the shortest projection, m = iV/8 = 32, 
the self-optimized state gives a projected energy which deviates by only 0.06% 
from the exact ground state energy, and for larger m the energy is exact within 
statistical errors. The other two trial states also lead to the correct energy 
but only for much larger m. The energy of the self-optimized trial state itself 
is shown in the right panel — its error is as small as that of the best variational 
amplitude- product state [33]. 
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Fig. 6. Long-distance spin correlations in the same runs as in Fig. 5. A different 
"random projected" state is picked for each m, which leads to an un-smooth curve. 
Fluctuations beyond the small statistical errors in the right panel reflect differences 
in how closely the rather short self-optimization runs have approached Pq — P m 



Fig. 6 shows the long-distance spin correlation calculated in the same 
runs. Again, the self-optimized state delivers superior results, although here 
the convergence is not as fast as for the energy. The error of the spin correla- 
tion in the trial state (right panel) is about 4% for large m, which is twice the 
error in the best variational state [33]. Thus, the self-optimized state is not 
identical to the best variational state, and an even faster convergence could 
be achieved by using the fully optimized variational state. However, the vari- 
ational calculation is much more time consuming than the self-optimization. 

We can go beyond the amplitude-product state by taking into account 
bond correlations. We are currently exploring this with both variational and 
self-optimized states. 



4 Triplet excitations 



A unique advantage of the VB basis is that an m z = triplet state can 
be projected simultaneously with the singlet, with essentially no additional 
overhead. Any triplet can be expanded in VB states where one of the bonds 
corresponds to a triplet [18]; ~ * where 

M = (I TiW + UiTi>)/>/2. (22) 
Formally, such a triplet can be created by acting on a singlet with Sf — 5|; 

(S?-S*)(i,j) = [i,j}. (23) 
To create a triplet |r(q)) with some momentum q, we can apply 

S* = -^Ve iq - r >S| (24) 

q Vn ^ 2 

3 
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(25) 



The amplitude-product state (18) for a periodic-boundary system has q = 
if the number of bonds, N/2, is even, whereas for odd N/2 it has q = (n, ir). 
This simply follows from the fact that sublattice A — > B and B —* A when 
translating by one lattice constant, whence each singlet (3) acquires a minus 
sign. We typically work with systems with even N/2 (e.g., Lx L lattices with 
even L) and so we will here consider q = singlets. 

At the antiferromagnetic wave- vector, Q = (ir,n), Sq acting on an ar- 
bitrary VB basis state \V) can be written as a sum of N/2 terms of the 
form (Sf — Sj)\V), with i, j corresponding to the sites connected by bonds. 
Operating on a q = singlet thus gives 



where the unspecified coefficients fk, e.g., the amplitude products in (18), 
have translational invariance built in. Thus, for a triplet with q = Q the 
wave function phases are buried in the definition of the singlets. Often, the 
lowest excitation of a Heisenberg system is a q = Q triplet, which we thus 
can sample without any difficulties with signs or phases. We consider this 
case first, before turning to triplets with arbitrary momentum. 

There are two possible actions of a singlet projector on a triplet bond; 



i.e., a diagonal operation on a triplet bond destroys the state whereas an 
off-diagonal operation on one triplet and one singlet bond creates a singlet 
at the sites on which the operator acts and moves the triplet to the other 
two sites involved. Importantly, the matrix element remains the same as in 
the off-diagonal operation on two singlet bonds. The weight of a triplet path 
is therefore the same as the corresponding singlet path (11), except that the 
triplet dies (giving zero weight) if an operator in P r acts diagonally on the 
triplet bond. We can thus measure triplet properties using paths generated in 
a singlet simulation, by considering only those triplet paths that do survive 
the propagation. In the trial state we have N/2 possible locations of the 
triplet. All of them can be attempted collectively in a single propagation, by 
keeping counters t(i) for the number of surviving states in which the triplet is 
connected to site i (with i on sublattice A). Initially t(i) = 1 for all i. During 
the propagation, for each diagonal operation (27) t(i) — > t(i) — 1, and for 
each off-diagonal operation (28) t(i) — ► t(i) — 1, t(l) — > t(l) + 1. Eventually, 



•SqK )) = SQ^ j fk\{i\,jl)k{i2,h)k ■ ■ ■ {iN/2,3N/2)k) 




(26) 



Hi 



M---> = o, 

[i,k]---{l,j)---) = \\---{i,3)---[l,k]---) 



(27) 
(28) 
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as m — > oo, all triplets die; t(i) = for all i, but typically there are enough 
survivors left at large enough m to compute converged triplet properties. 

An added advantage of calculating singlet and triplet properties in the 
same run is that there are error cancellations which in some cases can increase 
the statistical precision of differences, e.g., the singlet-triplet gap [36], 

A = E t (7t,tt)-Es(0,0), (29) 

by up to orders of magnitude relative to two independent calculations. The 
triplet energy Et(tt, 7r) can be estimated using (15), taking into account that 
a diagonal operation on a triplet bond gives zero, i.e., n& — > n& — n t , where, 
for surviving triplet configurations, nt = 0, 1 is the number of triplet bonds 
of length 1. Other triplet properties have been discussed in Rcf. [40]. 

We now discuss calculations with triplets of arbitrary momentum q. The 
energy can be evaluated according to 

_ E>(Q)l^ q gp^KQ)) 
E{ti) - E>(0)|SV^K0)> ' (30) 

We want to evaluate this expression using the sampled q = singlet bond 
configurations, and so we rewrite it as 

£Xo)|s* q ffP r s*Ho)> /£>(o)|s* q p r s*Ko)) 



{ ' EMWPrWm ~ \ £>(o)|p>(o)> J ■ [ ' 

The two factors can be evaluated based on sampling the propagations P r and 
the amplitude-product state |cr(0)) [which we have not explicitly written as 
a sum of bond configurations in (30) and (31)]. In the singlet energy (12), we 
could pick the Neel state for (if - ! and then obtained the very simple expression 
(15). Now we must consider the overlap with a momentum q triplet state. 
We use (cr(0)|S'i q , but in (31) have rewritten E(q) so that only the overlap 
with (cr(0)| has to be considered for the sampling weight. Phases arising from 
Sq only appear in the measurements, but in the end we have to evaluate 
the ratio of the two quantities in (31), which can be challenging in practice. 
However, close to q = (tt,it) and (0,0) [q ^ (0,0)] we find that it can be 
done; in some cases the method works even far away from these momenta. 
We define the dispersion relative to the gap (29) at (ir, w); 

w(q) - E T (q) - A. (32) 

Fig. 7 shows results for A and w(qi), where qi is the momentum closest to but 
not equal to (n, 7r); qi = (tt — 2tt/L, it). We show the convergence as a function 
of N/m for 4 x 4 and 16 x 16 lattices, comparing with exact diagonalization 
results in the former case. From (5) one would expect the convergence to 
be asymptotically exponential, which is seen clearly for L = 4. For L = 16, 
the three largest-m points are equal within statistical errors, suggesting that 
these results are also close to converged. Using w(qi) = 0.62 for the L = 16 
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system gives the spin- wave velocity c = 1.58, in very close agreement with 
the known value [43]. 

Another useful quantity accessible with the VB projector is the matrix 
element (r(q)|Sq|<7(0)), the square of which gives the single-magnon weight 
in the dynamic structure factor (which is experimentally measurable using 
neutron scattering). It can be calculated in a way similar to E(q), and we have 
done so successfully for q close to (ir, n). Results will be presented elsewhere. 

Finally, we also note that the triplet bond-length distribution gives a di- 
rect, albeit basis dependent, window into the "spinon" aspects [44] of the 
excitations. Spinon deconfinement should be manifested as a delocalized dis- 
tribution function, whereas two spinons bound into a magnon or "triplon" 
should be reflected in a well-defined peak in the distribution function. We are 
currently exploring this. 
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